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ABSTRACT 

We present a suite of cosmological radiation-hydrodynamical simulations of the assem¬ 
bly of galaxies driving the reionization of the intergalactic medium (IGM) at z > 6. 
The simulations account for the hydrodynamical feedback from photoionization heat¬ 
ing and the explosion of massive stars as supernovae (SNe). Our reference simulation, 
which was carried out in a box of size 25 h~^ comoving Mpc using 2 x 512^ particles, 
produces a reasonable reionization history and matches the observed UV luminosity 
function of galaxies. Simulations with different box sizes and resolutions are used to 
investigate numerical convergence, and simulations in which either SNe or photoion¬ 
ization heating or both are turned off, are used to investigate the role of feedback from 
star formation. Ionizing radiation is treated using accurate radiative transfer at the 
high spatially adaptive resolution at which the hydrodynamics is carried out. SN feed¬ 
back strongly reduces the star formation rates (SFRs) over nearly the full mass range 
of simulated galaxies and is required to yield SFRs in agreement with observations. 
Photoheating helps to suppress star formation in low-mass galaxies, but its impact on 
the cosmic SFR is small. Because the effect of photoheating is masked by the strong 
SN feedback, it does not imprint a signature on the UV galaxy luminosity function, 
although we note that our resolution is insufficient to model star-forming minihaloes 
cooling through molecular hydrogen transitions. Photoheating does provide a strong 
positive feedback on reionization because it smooths density fluctuations in the IGM, 
which lowers the IGM recombination rate substantially. Our simulations demonstrate 
a tight non-linear coupling of galaxy formation and reionization, motivating the need 
for the accurate and simultaneous inclusion of photoheating and SN feedback in models 
of the early Universe. 
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1 INTRODUCTION 

The first sources of ionizing radiation in our universe 
transformed the cold and neutral cosmic gas that was 
present shortly after the Big Bang into the hot and ion¬ 
ized plasma that we observe in the intergalactic medium 
(IGM) today. Research into this transformation, which 
took place in the first billion years, during the Epoch of 
Reionization, is crucial for understanding the formation 
and evolution of galaxies, includin g our own galaxy, the 


and evolution ot galaxies, includin g our own galaxy, tne 
Milky Way (for reviews see, e.g ., iBarkana Sz LoebI |200ll: 
Morales fc Wvith^ I, UUotti gORjP ^omm fc Yoshida 


20 111 : Zaroubi 2013l ~ Nataraian fc Yoshidal 12014 '). A num¬ 


ber of new telescopes are currently being developed or 
are already underway to unravel the astrophysics of these 
early times, including the Low Frequency Array (LOFAR; 
e.g., IZaroubi et al.l l2012li. the Mu rchison Widefield Ar¬ 


ray (MWA; e.g.. iLidz et_ajj _ 200^, t he Sq uare Kilome- 


Webb Space Telescope (J WST; e.g.. 


ter Array (SKA; e.g.. iMellema et ah 2013 ), the Jame s 


Zackrisson et al.l 
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and many others (e.g . . iFurlanetto. Oh. fc Briggsl 200d; 
iFan. Cariili. fc Keating l2006l : IPritchard Sz Loed 1201^ . 
These investments at the observational frontier have spurred 
the construction of models of the early Universe to guide the 
design of observational campaigns and to help interpret the 
observations once the data arrives. 

Most of the current models are built on the 
observationally supported notion that reionization 
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Choudhurv & Ferraral 2007: Volonteri & Gnedin 20091: Loeb 

2009 

• 

R.obertson et al. 201fll: Raicevic. Theuns. & Lacev 

2011 

Fontanot. Cristiani. & Vanzellal 2012|). A natural 


implication of these models is a strong coupling between 
reionization and galax y formation caused by ionizing 
radiative feedback (e.g.. iQardi &i Ferrarall2005l l. Photoion¬ 
ization heating by star-forming galaxies typically raises the 
temperature in the reionized gas to a few times K. 
The associated increase in gas pressure increases the Jeans 
scale, i.e., the scale below which pressure prevents the 
collapse of ga s into gravitationally bound objects (e.g., 
lGnedinll2000bt ) . This boils gas out of low-mass dark matter 
(DM) haloes and im pedes the accretion of gas from the 
IGM onto them fe.g..lR eeslll98 6l:ISha,r )irOjGirouXj_&_^abiJ 


UjiVi onto tnem (e.g., IK gegLj^ool;lM^BinhJjirQUx^_fe^DM)U: 

1994; Barkana fc Loe^ 199^ Gnedin l2000b : Diikstra et al 


2004 1 _ FinlatorjDavejOze] 2011 : IPetkova fc Springell 201ll : 
Noh fc McQuinnll2014l ~).~Pliotoheating therefore leads to a 
reduction in the fuel from which stars form, which lowers 
the star formation rates (SFRs) and the ionizing emissivi- 
ties of low-mass galaxies. This makes it more difficult for 
galaxies to reionize the IGM and, therefore, photoheating 
provides a negative feedback on reionization. 

Photoheating also provides a positive feedback on reion¬ 
ization. The increase in the gas pressure and Jeans mass 
smooths density fluctuations in the IGM. This decreases 
the gas clumping factor, Cigm = (n^)/?!^, where the an¬ 
gular brackets indicate the average in t he IGM, and hn is 


Pawlik. Schave. & van Scheroenzeel 


2009 

: IFinlator et al.l 

2 OI 2 I: Shull et all 20121: ISo et al. 


2014 

). Because the 

rate at which the IGM recombines is 

proportional to 


the gas clumping factor (e.g., iMadau. Haardt. fc ReesI 
Il999h . photoheating decreases the IGM recombination 
rate. Fewer ionizing photons are then needed to keep 
the gas ionized, which facilitates reionization (e.g. , 

Pawlikj_SchaWj_^van_Scher£enze^ 20091: Ishull et al1l20ld : 


feed back on the formation of structure and the reionization 
of the IGM makes modelling reionization a challenging task. 

Cosmological simulations are generally considered to 
be among the most powerful techniques to investigate 
the impact of s tellar feedback on rei o nization and g alaxy 
formation (e.g., iTrac fc GnedinI l201ll : iFinlatoJ 12012 1 , but 
are also computationally highly demanding. A princi¬ 
pal computational challenge is the required large dy¬ 
namic range. The characteristic size of individual ion- 
ized regions is of order 10 comoving Mpc (cMpc; e.g., 
iFurlanetto. McQuinn. fc Hernauisd l2006ll . The redshift at 
which these regions percolate to reionize the Universe shows 
spatial fluctuations due to the modulation of galaxy forma¬ 
tion by structure formati on modes on still larger scales, of 
order 100 h~^ cMpc le.g.. lBarkana fc Loeb|[^04l : llliev et al.l 
l2014h . The concurrent need to resolve the Jeans mass in 
the IGM as well as the population of atomically cooling 
galaxies with masses above 10® M©, the prime candi¬ 
dates for driving reionization, extends the sp ectrum of rele¬ 
vant scales down to a few kp c, and less le.g.. iGnedin fc FanI 
I2OO6I : iBolton fc Beckerj|20od ~l . 

A second computational challenge is the large ex¬ 
pense of the radiative transfer (RT) of ionizing photons. 
A straightforward and accurate treatment of the RT con¬ 
sists of tracing photons from each ionizing source along 
straight rays, one sour ce at a time (for an overv i ew of 
RT methods s ee, e. g., Illiev et ^ I I 2 OO 6 I : Illiev et all l2009l : 
iTrac fc GnedinI I 2 OIil l. The computational cost of such ap¬ 
proaches therefore increases linearly with the number of 
ionizing sources. This renders simulations of cosmological 
volumes that contain large numbers of sources computa¬ 
tionally demanding. The issue becomes particular severe 
if the ionizing radiation f rom recombining ions needs to 
be treated accurately (e.g., Mjaseni^_Fernaxa ^_fcCiardill2003l : 
iHasegawa fc Umemura|[2010l : lRaicevic et al.li2014l ~). The high 
cost of ray-based RT methods has triggered the development 
of moment-based RT methods that ar e computationally less 


Finlator et al.| |2012|; |Sobacchi fc Mesingeij |2014|) . This pos- expensive, but also less accurate (e.g., Gnedhi&_^be| 200^ 


itive feedback from photoheating will be especially strong 
if X-ray sources preheat the gas before it is reion- 


ized (e.g.. Madau et a 

.1 2 OO 4 I: R cotti. Ostriker. & GnedinI 

200Sl: iHaimanl 1201 ll: 

Jeon et al. 2014al: Xu et al.l 20141: 

Knevitt et al.ll2014h. 

n photoheating may be accompa- 
the explosion of stars as supernovae 
^Il98fil: IWvithe & Loebll2013lL SNe 
s star formation in low-mass galax- 
winds, which makes it more diffi- 

The feedback froi 
nied by feedback from 
fSNe: e.g.. iDekel & Sil 
are thought to suppres 
ies by expelling gas in 


Aubert fc Tevssierll2010l:lFinlator et al.l 20091 : Rosdahl et al.l 

2 OI 3 I : Norman et al. 2013l h 


To simplify the problem, cosmological hydrodynamical 
simulations of galaxy formation often assume the gas to be 
in ionization equi librium with a spatially uniform ionizini 


background (e.g., Miralda^EscudejIIaehneltj_&R^^ 200 


Pawlik. Schave. fc van Scherpenzeell 20091 : Duffy et al.l 


201J). Such simulations enable the investigation of the 


cult for galaxies to reionize the IGM. However, the ejection 
of gas may open additional low-density channels through 
which ionizing photons may escape the ISM more easily, 
thus making i t easi e r for galaxies to re i onize the IGM (e.g.. 


response of the gas in galaxies and the IGM to heating by 
ionizing photons at relatively low computational cost, but 
must assume the redshift at which the ionizing background 
is instantaneously turned on. The main drawbacks of such 
simulations are that they do not yield insights in the 
timing and morphology of reionization, and that fluctu- 


Yaiima et al.l 2003: Wise Sz Genll2009l: IPaardekooper et al.l at ions in the intensity due to local sources (e.g., Gnedin 


20 111 : Kimm fc Cen 20141 '). The suppression of star for¬ 
mation by SNe may be amplified in a nonlinear man- 
by the feedback from photoheating (|PawlilL_fc_Schay^ 


ner by tne teed pack trom pnotoneatmg l|Jr 'awlik: t z bcnavel 
I2OO9I : iFinlator. Dave. Ozell I2OIII : iHopkins et al.l 2013l l . Lo¬ 
cally, SNe may also increase the rates at which stars form 
as interstellar gas is comp ressed to star-forming densitie s 
in colliding SN shocks (e.g.. lGeen. Slvz. fc Devri^dtll2013l) . 
The wide range of possibilities by which star formation may 


lYaiima. Choi, fc Nagamind l2012l : iRahmati et al 


2013bl ') and the ability of the gas to throw shadows and 


self-sh i eld (e.g., Abel. Wise, fc BrvanI l2007l : IRahmati et al.l 
l2013al: lAltav et all I 2 OI 3 II are ignored. Sometimes, self¬ 
shielding is approximately accounted for by setting the 
intensity of the ionizing background to zero at densities 
above which the gas is expe cted to be neutral in the 
absence of local sources (e.g., iNagamine. Choi, fc Yaiimal 
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2010 I: IPaardekooper, Khochfar. fc Dalla Vecchial l2013l : 

Vogelsberger et al. 2013l l. 


Other simulations focus on the accurate treat¬ 
ment of the RT of ionizing photons at the expense 
of a detailed treatment of the feedback from radia- 


tive heatine (e.g.. Nakamoto. Umemura. & Susal 

2 OO 1 I 

Razoumov et al.l l2002l: ICiard 

. Stoehr. & Whitel 

20031 

Iliev et al.l 20061: McOuinn et al.l 

2007: Finlator et al.l 

200d 

Aubert fc Tevssied 120101: iBaek. Ferrara, fc SemelinI l2012f) 


In this type of simulations, the computational complexity 
is reduced by computing the RT and the evolution of 
structure separately, ignoring the impact of photoionization 
heating on the dynamics and distribution of the gas. Often, 
to simplify the problem further, the dynamics of the gas 
is ignored altogether and assumed to trace that of the 
DM ('e.g.. Illiev et al.ll2006l i. A semi-analytical model then 
specifies the location of the galaxies and the rates at which 
these galaxies form stars. Sometimes, the semi-analytical 
model is calibrated to approximately account for the 
radiation-hydrodynamical respons e to photoheating (e.g., 
ICiardi et ai.ll20o4 Illiev et ^120071 ). 

The most direct but also computationally most expen¬ 
sive way to investigate the coupling of reionization and 
galaxy formation, is to carry out three-dimensional cos¬ 
mological radiation-hydrodynamical simulations. Several of 
th ese simulations now exist, inclu ding the pionee ring works 
bv iGnedin fc Ostrikeil iIIQQTII and iGnedinI (I 2 OOOII . However, 
even after more than a decade of intense research, the field 
is still in its infancy. Radiation-hydrodynamical simulations 
have only very recently begun to approach the cosmological 
scales an d the high resolution needed to capture the relevant 
phys ics jNorman et al.l l2013l : iGnedinI l2014l : iKimm fc GerJ 
I2OI4I) . 

Most earlier radiation-hydrodynamical simulations 
of reionization were restricted to relatively small vol- 


cal mean (e.g.. ICnedinl l2000b 

: iRicotti. Gnedin, fc Shulll 

20081: Petkova & SDringelll201ll: 

Finlator. Dave. Ozell 201ll: 

Hasegawa fc SemelirJl2013|). Onlv some simulations afforded 


ray-based RT methods that ensure the accurate propaga- 
tion of ionization fronts in inhomogeneous density fields 
dHasegawa fc SemelinI [^13ll . The majority of simulations 
employed the computationally less expensive but also nu¬ 
meric ally diffusive momeri t-base d RT methods (e.^., Gnedd 
2 OOOI: Petkova fc S pringell 201ll: Finlator. Davi. Ozelll201ll : 


Norman et al. l l2013l : ISo et al.ll2014ll . In some cases, the com¬ 


putational expense needed to be reduced further by lower¬ 
ing the spatial resolution at which the RT is carried out, 
thus requiring a sub-resolution model for the propagation 
and co nsumption of photons at h ydrodynamically resolved 
scales dFinlator. Dave. Ozell[201ll '). 

In this work we present radiation-hydrodynamical sim¬ 
ulations of reionization in cosmological volumes at high spa¬ 
tially adaptive resoluti on, implementing accur ate RT using 
the TRAPHIC code JPawlik fc Schavd l2008lb TRAPHIG 
solves the time-dependent RT equation by tracing photon 
packets at the speed of light and in a photon-conserving 
manner through the simulation box. TRAPHIC thus en¬ 
ables an accurate treatment of shadows and self-shielding. 
The photon packets are transported directly on the spa¬ 
tially adaptive, unstructured grid defined by the SPH par¬ 
ticles, thus exploiting the full dynamic range of the hydro¬ 


dynamic simulation. The large number of ionizing sources 
typical of these simulations does not pose a problem be¬ 
cause TRAPHIC employs a photon packet merging tech¬ 
nique that renders the computational cost of the RT in¬ 
dependent of the number of sources. TRAPHIC is imple¬ 
mented i n the galaxy f ormation code GADGET (last de¬ 
scribed in ISDring^l2005l l . and this enables us to account for 
the radiation-hydrodynamical feedback from star formation 
and reionization. 

We focus our investigations on the roles of radiative 
heating (which we also refer to as photoheating) and SNe 
in shaping the SFR and reionization history in the phys¬ 
ically motivated and observationally supported models of 
the early Universe that our simulations provide. It is cur¬ 
rently not feasible to provide an ab initio description of 
the internal structure of galaxies in cosmological volumes, 
even in substantially less expensive simulations in which 
the detailed transfer of ionizing photons is ignored (e.g., 
Trac fc GnedirJli^OllI : IScannapieco et aDl2012l : ISchave et al'l 
2015h . Simulations like ours therefore rely on sub-resolution 
models of the interstellar medium (ISM). This introduces 
free parameters that require careful calibration, e.g., using 
zoomed simulations that achieve a higher resolution by fo¬ 
cusing on a smaller computational volume, or observations. 
The two main parameters relevant to simulations of galaxy 
formation during reionization are the fraction /sn of the en¬ 
ergy released in SNe and other unresolved stellar evolution 
processes that is available to heat the ISM, and the fraction 
js^bres ionizing photons that escapes the unresolved 

ISM and are available to reionize the gas. 

This paper is organized as follows. In Section [5] we de¬ 
scribe our numerical techniques and the set of simulations 
we have carried out. In Section[3]we describe our results. We 
start in Section ED with an overview of the simulations and 
a comparison with observational constraints. This is contin¬ 
ued in Section [3. 2 1 and extended by an analysis of the impact 
of stellar feedback on galaxy formation and reionization. In 
Section E we discuss the sensitivity of our results to vari¬ 
ations in physical parameters as well as some of the main 
limitations inherent to our numerical approach. In Section [5] 
we conclude with a brief summary. Lengths are expressed in 
physical units, unless noted otherwise. 


2 SIMULATIONS 

2.1 Gravity and hydrodynamics 

We use a modified version of the N-body/TreePM S moothed 
Particle Hydrodynamics (SPH) code GADGET (ISpringell 
I 2 OO 5 I: our specific ve rsion is derived from that discussed 
in ISchave et al.l l20ldj to perform a suite of cosmological 
radiation-hydrodynamical simulations of galaxy formation 
down to redshift 2 = 6 (see Table [T]). 

The simulations are initialized at redshift 2 = 127. Ini¬ 
tial particle positions and velocities are obtained by applying 
the Zeldovich approximation (Zeldovich 1970) to p articles 
arran ged along a uniform grid of glass-like structure (IWhitd 
Il996ll . We use a transfer function for matter p erturbations 
generated with GAME llLewis fc Bridle! l200d l and apply 
it to describe perturbations in both the dark matter and 
the gaseous components. We adopt the ACDM cosmological 



















































































































4 Pawlik et al. 


Table 1. Simulation parameters: simulation name; comoving size of the simulation box, -Lbox! number of DM particles, mass of 

DM particles, moM; mass of gas particles, mgas; comoving gravitational softening scale, tgoft; purpose. The number of SPH particles 
initially equals A^dm (it decreases during the simulation due to star formation). All simulations are carried out down to redshift z = 6, 
except for simulation L12N512, which was stopped at 2 = 6.8 because of the large computational expense. 


simulation 

-^box 

[ h~^ cMpc] 


ruDM 

[Mq] 

m.gas 

[Mq] 

^soft 

[ h~^ ckpc] 

purpose 

L25N512 

25.0 

512 

1.00 X lO *" 

2.04 X 10® 

1.95 

reference 

L12N256 

12.5 

256 

1.00 X 10 ’^ 

2.04 X 10® 

1.95 

small box 

L12N512 

12.5 

512 

1.25 X 10® 

2.55 X 10® 

0.98 

high resolution 

L50N512 

50.0 

512 

8.02 X 10 ’^ 

1.63 X 10 ’' 

3.91 

large box 

L25N256 

25.0 

256 

8.02 X 10 ’^ 

1.63 X 10 ’' 

3.91 

low resolution 


model with parameters = 0.265, fib = 0. 0448 and 12 a = 
0.735 , Us = 0.963, as = 0.801, and h = 0.71 llKomatsu et al.l 
[2 oT 3). These parameters are consistent with the most re¬ 
cent constraints from observations of t he Cosmic Microwave 
Background by the Planck satellite Jpianck CollaborationI 
l2014a : IPlanck Collaboration et al.ll2015^ . 

Our reference simulation utilizes 2 x 512^ DM and gas 
particles in a box of size 25 h~^ Mpc. This implies a DM 
particle mass of mum = 0.7 x 10^ h~^ Mq « 10^ M© and a 
gas particle mass of rugas = 1.4 x 10® h~^ Mq « 2 x 10® MfT) . 
Atomically cooling haloes with mass (Eq. 3.12 in lLoe 


Mvir ^ 10 Me 


104 K 


3/2 


. 6 / 


3/^ fl + z 

TiT 


-3/2 


( 1 ) 


where /r is the mean atomic weight, are therefore resolved by 
at least Mvir(Dm — ~ 10[(1 -I- 2)/10]~®/^ DM 

particles. 

We apply a comoving Plummer-equivalent gravitational 
softening length, Cgoft, equal to 1/25 of the mean DM par¬ 
ticle separation to both DM and baryonic particles, i.e., 
Csoft « 1.95 h~^ ckpc (I//25.0 h~^ cMpc) (A7j!,{^/512)“^, 
where A^dm is the number of DM particles. SPH quanti¬ 
ties are estimated by averaging inside a sphere containing 
fVngb = 48 neighboring gas particles and adopting the en- 

S conserving formulation of SPH (ISoringel fc HernauistI 
. The SPH kernel, i.e., the radius of this sphere, is pre¬ 
vented from falling below 10“^esoft. 


2.2 Chemistry and cooling 


While our simulations account for the thermal feedback from 
the explosion of massive stars as core-collapse SNe (see Sec¬ 
tion EUl, for simplicity, we do not follow the ejection and 
transport of metals. We therefore follow the non-equilibrium 
chemistry and radiative cooling of the gas, assuming that 
the gas is of primordial composition. Molecule formation 
is assumed to be suppressed by the Lyman-Werner back¬ 
ground expected to pervade the u niverse soon after the 
form a tion of the first stars (e.g.. iHaiman. Rees, fc Lo^ 
I 997 I: Ricotti. Gnedin, fc Shulll 2001 : Wise fc Abd 20051 : 


Greif fc BrommI 120061 : I Aim et al.l 120091 ). a process not re¬ 
solved in our simulations. Our simulations thus ignore that 
molecular hydrogen may form in self-sh ielded regions (e.g., 
IWolcott-Green. Haiman. fc BrvarJl201ll ). 

We consider all relevant ato mic radiative cooling pro - 
cesses, using the rate coefficients in IPawlik fc Schavd (l201ll) : 
cooling by collisional ionization, collisional excitation of 


atomic lines, the emission of free-free and recombination ra¬ 
diation, and Compton cooling by the CMB. The species frac¬ 
tions and the gas temperature are advanced in time using 
the explicit subcycling t ime integration scheme presented in 
IPawlik fc Schava (l201ll ). Once stars form and emit radia¬ 
tion, the chemical and thermal evolution of the gas is also 
affected by photoionization and photoheating, and by the 
injection of energy by SNe, as we describe below. 

Solving for the chemistry and cooling of the gas is by far 
the most expensive operation in our simulations. To speed 
up the calculations and unless otherwise noted, we set X = 
1, i.e., we neglect helium. In Section l4.ll we will describe 
a preliminary comparison with a simulation in which we 
follow the radiative heating and cooling by helium assuming 
X = 0.75. 


Our simulations do not have sufficient resolution and 
do not capture the physics required to describe the multi¬ 
ple phases of the dense gas in galaxies that are expected 
to develop a bove densities nn ~ 10“^ — 10~^ cm~® (e.g., 
ISchaveir2004h as observed in the local universe. We therefore 
employ a subgrid model to describe the thermodynamics of 
the dense gas, imposing a pressure floor P oc on gas 
with densities nn > nfj, where =0.1 cm“®, normalized 
to P/k = 10® cm“® K at the critical density We use 
7 efr = 4/3 for which both the Jeans mass and the ratio of 
the Jeans length and the SPH kernel of gas at the pressure 
floor are independent of the density, thus preventing spu- 
rious fragmentation due to a lack of numerical res olution 
(iBate fc Burkertlll997l : ISchave fc Dalla Vecchiall2008l) . Note 
that, unlike in simulations that impose a polytropic equa¬ 
tion of state, in our simulations, gas with densities above the 
critical density can have pr essure above the pressur e floor . 
Our implementation follows Sch ave fc Dalla Vecchial (l2008l) 


and IPalla Vecchia fc Schave ( 20121 ). to which we refer the 
reader for further details. In our simulations adopting a hy¬ 
drogen mass fraction X = 1, we use a critical density higher 
than riH by a factor 1/0.75 to facilitate comparisons with 
simulations in which X = 0.75. 


2.3 Star formation 

We employ the sta r form ation implementation of 
ISchave fc Dalla Vecchial ll2008l) . to which we refer the 
reader for details. Briefly, gas with densities nn > 
and temperatures T < T* is allowed to form stars using 
a pressure-dependent rate that reproduces the Kennicutt- 
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Schm idt relation observed in the local universe dKennicutd 
Il998h in idealized simulations of isolated disk galaxies, 

E* = 1.515 X 10-^ M 0 yr kpc'^ Me^pc-^ ’ 

where E* is the SFR surface density and Egas the gas surface 
density. The last equation has already been renormalized by 
a factor o f 1/1. 65 to account for the fact that it assumes the 
ISalpeterl lll955ll IMF whereas we are using a Chabrier IMF. 
This co nversion factor between SF Rs has been computed us¬ 
ing the lBruzual fc Charlotl (l2003l i population synthesis code 
for model galaxies of age > 10 ^ yr forming stars at a con¬ 
stant rate and is insensitive to metallicity. 

We set Uh = and use T* = max(10® K, lO'^ ^Tfloor), 
where Tfioor is the density-dependent temperature implied 
by the pressure floor below which particles are prevented 
from cooling. Moreover, to avoid numerical artefacts at high 
redshifts when the cosmological mean gas density is high, we 
require a minimum overdensity of 57.7 for gas to be allowed 
to form stars. 

The star formation law is interpreted stochastically, and 
the probability that a star-forming gas particle is turned into 
a star particle in a time interval At is given by min( At/r*, 1), 
where r* = Egas/E* is the gas consumption time. Gas par¬ 
ticles are converted to star particles assuming a conversion 
efficiency of 100 %, i.e., the masses of the star particles are 
identical to those of the gas particles from which they are 
formed. At higher resolution, star formation and the associ¬ 
ated feedback may be more extended in time and in space. 
This may affect the impact of star formation on the gas and 
therefore alter, e.g., the resolved escape fraction of ionizing 
photons. 


2.4 Ionizing luminosities 

We interpret the star particles as simple stellar populations, 
i.e., coeval stellar clusters that are characterized by an initial 
mass function (IMF), metallicity, and age. 

We compute the time-dependent hydrogen ionizing lu¬ 
minosities of these star fo rmation bursts using t he popu¬ 
lation synthesis models of ISchaererj (l2003ll . Since l^haereil 
ll2003l) does not tabulate ionizing luminosities for the 
Chabrier IMF, which is the IMF used by our star f ormation 
recipe , we adopt the ionizing luminosities of the ISchaereil 
(l2003tl models with Salpeter initial mass function in the 
range 1 — 100 Mq and metallicity Z = 0.02 = Zq. This 
yields ionizing luminosities close to those of stellar popula¬ 
tions drawn from a Chabrier IMF in the range 0.1 — 100 Mq 
with subsolar metallicity Z — 10“® = 0.05 Zq and is consis¬ 
tent with the IMF used in our star formation rec ipe and ap¬ 
propriate for stellar populations at 2 > 6 . fe.g.. iMaio et al.l 
120101 : Iwise et al.ll20l4[Jeon et al.]l2014bf l. 

We multiply the luminosities of each star particle by the 
sub-resolution escape fraction to model the removal 

of photons due to absorption in the unresolved ISM. Some 
previous works that carried out RT simulations of reioniza¬ 
tion have calibrated the ISM escape fraction to yield reion¬ 
ization histories i n agreement with curren t observational 
constraints (e.g., lAuber t fc Tevssiej l2010l : iFinlator et al.l 
l2012l : [Ciardi et ^l2012ll . However, as we will discuss below, 
the interpretation of these constraints is subject to substan¬ 


tial systematic uncertainties and limited by our incomplete 
understanding of the physical processes at play. We therefore 
prefer to set = 1 > a-nd use our simulations to explore 

the dependence of the reionization history on box size and 
resolution. In Section 1441 we will briefly discuss the impact 
of variations away from = 1 on our results. 


2.5 Ionizing radiative transfer 

We use the RT code TRAPHIC llPawlik fc Schav3 l2008l . 
I 2 OIIII to transport ionizing photons. TRAPHIC solves the 
time-dependent RT equation in SPH simulations by trac¬ 
ing photon packets emitted by source particles through the 
simulation box in a photon-conserving manner. The photon 
packets are transported directly on the spatially adaptive set 
of SPH particles and hence the RT exploits the full dynamic 
range of the hydrodynamical simulations. A directed trans¬ 
port of the photon packets radially away from the sources is 
accomplished despite the irregular distribution of SPH par¬ 
ticles by guiding the photon packets inside cones. A photon 
packet merging technique renders the computational cost 
of the RT independent of the number of radiation sources. 
In the following, we provide a brief overview of TRAPHIC 
in order to motivate the meaning of the numerical parame¬ 
ters of the RT s pecified below. The read er is referred to the 
description s in jPawlik fc Schavd (| 2008l'). Pawlik fc Schavel 
ll 201 lll and IPawlik. Milosavlievic, fc BrommI ( 2013l ~) for fur¬ 
ther discussion. 


2.5.1 Method 

The transport of radiation starts with the emission of pho¬ 
ton packets by star particles in N^c tessellating emission 
cones. The photons in each photon packet are shared among 
the iVngb < Angb neighboring SPH particles residing in the 
cones. In cones containing zero neighbors, an additional, so- 
called virtual particle is inserted to which the photon packet 
is assigned. Star particles emit photons using emission time 
steps Atem, in between which the orientation of the cones 
is randomly changed to increase the sampling of the volume 
with photons. The spectrum of the emitted radiation is dis¬ 
cretized using Nu frequency bins, and each photon packet 
carries photons from one of these bins. Therefore, the num¬ 
ber of photon packets emitted per emission time step is 
Nv X Abc- 

The newly emitted photons are assigned a propagation 
direction parallel to the central axis of the associated emis¬ 
sion cone, and, together with any other photon in the sim¬ 
ulation box, are then propagated to the downstream neigh¬ 
bors of the SPH particle at which they reside. A particle is 
a downstream neighbor if it is among the iVngb neighboring 
gas particles and resides in the regular transmission cone 
centred on the propagation direction and subtending a solid 
angle of dvr/Arc- The parameter Arc defines the angular 
resolution of the RT. If there is no downstream neighbor 
inside a transmission cone, a virtual particle is created to 
which the photons are then propagated. The transmission 
cones confine the propagation of photons to the solid angle 
in which they were originally emitted. The transport of pho¬ 
tons occurs at the user-specified speed c and is discretized 
using RT time steps Atr. 
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A given SPH particle may receive several photon pack¬ 
ets within the same RT time step At^. These photon packets 
are collected according to their propagation directions using 
a set of A^rc tessellating reception cones. Photon packets 
whose propagation directions fall in the same reception cone 
are merged and replaced by a single new photon packet. 
Each reception cone subtends a solid angle An/N- rc- Hence, 
the parameter Nrc determines the angular resolution of the 
merging. In the absence of virtual particles, the merging 
strictly limits the maximum number of photon packets in 
the simulation box to Nrc x Ni, x A^sph, where A^sph is 
the number of SPH particles, and renders the computation 
time independent of the number of sources. Note that the 
angular resolution at which photon packets are merged may 
be chosen independently of the angular resolution at which 
photon packets are transferred. 

Photons are absorbed as they propagate through the gas 
from SPH particles to their neighbors depending on the op¬ 
tical depth between the tw o neighboring particles, respect¬ 
ing photon conserv ation (lAbel. Norman. &: Madaul 1 19991 : 
iMellema et aLlEoOfil l. The absorption of photons within each 
frequency bin is treated in the grey approximation using 
photoionization cross section:{3 


((THl)i- = / dv ---CrHl(i^) X / 

j l/\ (17* h' _J jy 


du- 


hpu 


,(3) 


where Jv{v) is the spectrum, and u\ and Vh are the low 
and high energy limits of frequency bin v. The number 
of absorbed photons determines the photoionization rate 
P-vHi.i/ of HI in the given fr equency bin v defined by (e.g., 
lOsterbrock &: FerlandllioOfil l 


■p / \ j AnJv{v) , . 

r 7 Hi,i- = {om)v / dv — -. (4) 

The photoionization rate implies a photoheating rate given 
by — (£:Hi):/r 7 Hi,i^, where 




‘" 1 ' AtiJu{v) 

dv - 

hpv 


<Jm{v){hpv — hpvm) 


‘'i' j AnJ^{v) 

dv — - -crHi(i^) 

hpv 


(5) 


is the grey excess energy of frequency bin v, and hvRi = 

13.6 eV is the photoionization threshold energy of HI. The 
photoionization and photoheating rates are then passed to 
a chemistry solver to update the HI fraction and the tem¬ 
perature of the gas. 


2.5.2 Numerical parameters 

In the simulations presented in this work, we set Nngb = 
32, and choose an angular resolution of the transport of 
Ntc = 8 and of the merging of Arc = 8. Sources emit 
photons into Arc = 8 directions using emission time steps 
Atem = min(10“^ Myr, Atr). Photons are transported at a 
speed c = c, where c is the speed of light, using time steps 
of size Atr = min(10“^ Myr, Athydro), where Athydro is the 


^ While TRAPHIC can treat ionization of helium, in this work 
we set X = 1, as discussed in Section 12.21 


smallest GADGET particle time step. We use a single fre¬ 
quency bin with bounding energies located at [13.6, oo] eV. 
We use the fits to the frequency-dep endent photoionization 
cross sections bv IVerner et al.l (Il996ll and adopt a black body 
spectral shape of temperature Tbb = 5 x 10"^ K, consistent 
with the parameters of the adopted population synthesis 
models, to compute the grey photoionization and photoheat¬ 
ing rate associated with that bin. This gives (chi) = 2.93 x 
10“^® cm“^ and (ehi) = 3.65 eV. The RT assumes periodic 
boundary conditions, i.e., photons leaving the box through 
any given face are reinserted at the periodically opposing 
face. The de pendence of the RT on th e par ameters has been 
discus sed in lPawlik fc Scha^ ll2008ll and IPawlik fc Schavel 
ll201lll . to which we refer t he reader. Additional tests can be 
found in th e appendices of iPawli k. Milosavlievic. fc BrommI 
ll2013tl and iRahmati et aH i 2013a 1. For simplicity, we treat 
recombination radiation in the on-the-spot approximation, 
i.e., we compute recombination rates using Case B (e.g., 
iRaicevic et al.l[2014l : iTanaka et al.|[20lll . 


2.6 Energy injection by core-collapse supernovae 

The feedback from the explosion of stars as core-collapse 
SNe is modeled as an injection of thermal energy into the 
ISM surrounding the star particlefl Our numerical imple¬ 
mentation, which has been designed to control spurious ra¬ 
diative energy losses due to the limited numer ical resolution, 
is described in iDalla Vecchia fc Scha^ (l2012l l. Core-collapse 
SNe inject thermal energy after a delay of 30 Myr after 
the birth of the star particle, approximately correspond¬ 
ing to the maximum lifetime of stars that end their lives 
as core-collapse SNe. For each core-collapse SN that occurs, 
/sN X 10®^ erg of thermal energy is injected. The energy 
is distributed stochastically to an average subset of 1.34 of 
the 48 neighboring SPH particles, instantaneously increas¬ 
ing their gas temperature by AT — 10^ ® K and assuming 
a mean atomic weight p, = 0.6 appropriate for fully ionized 
gas. We set /sn = 1, which results in a good agreement of 
the star formation history in our reference simulation and 
observational constraints a.t z > 6 (see Figure [4]). 


2.7 Identification of galaxies 

We post-process our simulations to extract haloes using 
the friends-of-friends (FOF) halo finder, with linking length 
equal to l/5th of the mean DM inter-p article distance, built 
into the substruc ture finder Subfind (ISpringel et al.l 1200 ll : 
iDolag et al.ll200^ . For each FOF halo, we use Subfind to 
identify the particle for which the gravitational potential is 
minimum, and let its location mark the halo center. We ob¬ 
tain the virial radius as the radius of the sphere centered 
on the halo center within which the average matter density 


^ We refer to this thermal feedback as SN feedback for simplicity 
of presentation. However, our implementation is agnostic to the 
physical processes that inject the energy, and therefore may be 
interpreted as collectively describing all promptly acting feedback 
processes associated with the formation of stars that are not ex¬ 
plicitly treated otherwise. Such processes include st ellar winds as 
well as radiation pressure (for a discussion see, e.g. JSchave et ^ 
l2015ll . 
























































equals 200 times the redshift-dependent critical density of 
the universe. The total mass inside this sphere is the virial 
mass of the halo. The total SFR of the halo is the sum of 
the SFRs of the gas particles it contains. We will also com¬ 
pute the baryon fraction of the halo, which is the ratio of 
the total mass in gas and stars inside the virial radius and 
the virial mass. 


3 RESULTS 

In Section o we provide an overview of the simulations 
listed in Table [T] Thereafter, in Section [3.21 we discuss how 
stellar feedback from SNe and photoheating impact the for¬ 
mation and evolution of galaxies and the reionization of the 
gas. We will make use of additional simulations identical to 
those in Table [T] except that we have turned off SNe and/or 
photoheating to help isolate the impact of the feedback these 
processes provide. Most observational works on UV galaxies 
assume a Salpeter IMF to infer SFRs, while our simulations 
assume a Chabrier IMF. Where necessary, we have divided 
observationally inferred SFRs by a factor of 1.7 to account 
for the difference between the two IMFs (see Section [2.311 . 
Note that the high-resolution (high-res) simulation L12N512 
has been stopped at z = 6.8 due to its large computational 
expense. 


3.1 Overview of simulations 

Figure [T] shows images of the gas densities, temperatures 
and neutral and ionized hydrogen fractions in our reference 
simulation L25N512 at z « 9. By this redshift, reionization 
has already proceeded significantly, and individual reionized 
regions are just about to percolate the simulation box. 

Most of the gas inside the ionized regions is photo- 
heated to about 2 X 10^ K, the characteristic tempera¬ 
ture of stellar HII regions. Near galaxies, gas may be heated 
to much higher temperatures in gas-dynamical shocks that 
accompany the accretion of gas in massive haloes and the 
explosion of stars as SNe. Some of the gas in the ionized 
regions is sufficiently dense to shield from the ionizing ra¬ 
diation and remains neutral. The high resolution enabled 
by our spatially adaptive RT technique is key to tracking 
this gas phase. For comparison, a RT simulation on a uni¬ 
form grid with a number of grid cells equal to the number of 
SPH particles employed here would have limited our ability 
to resolve the structure of the ionized gas to scales above 
> 70 kpc comoving. This is larger than the typical scale of 
self-shielded (Lyman-limit) systems that are thoug ht to be 
among the main consumers o f ionizing photons (e.g ., ISchavel 
I 2 OOII : Ih^rlanetto &: OhilioOSi : ICnedin fc FanI 1200611 . Match¬ 
ing our spatially adaptive resolution using a uniform grid 
would require 512 x 25 = 12800 grid points per dimension, 
which is beyond current computational capabilities. 

Figure [ 2 ] compares a set of basic observables extracted 
from our simulations with current observational constraints. 
The SFR density of currently observationally accessible 
galaxies with SFRs ^ 0.3 Mq yr“^ (assuming a Chabrier 
IMF, see Figure |3|), is insensitive to changes in box size and 
resolution and in good agreement with observations. The 
comparison with the total SFR density suggests that cur¬ 
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rent observations uncover only a fraction of stars forming at 

There is a significant dependence of the total SFR on 
resolution. Because density fluctuations are increased at 
higher resolution, the reference simulation L25N512 yields 
initially a higher total SFR than the low-resolution (low- 
res) simulation L25N256, and a lower total SFR than the 
high-res simulation L12N512. On the other hand, near the 
final simulation redshift at 2 : = 6, the total SFR is slightly 
smaller in the reference simulation than in the low-res sim¬ 
ulation, and slightly larger than in the high-res simulation. 
This inversion in the dependence on resolution is caused by 
the stronger suppression of star formation by stellar feed¬ 
back at higher resolution. 

The top right panel of Figure [2] shows that the reion¬ 
ization history, i.e., the evolution of the neutral hydrogen 
fraction, is insensitive to the size of the simulation box 
L > 12.5 h~^ Mpc. This rapid convergence may partly re¬ 
sult from our use of identical amplitudes and phases that 
characterize the initial Gaussian random density fluctua¬ 
tions from which our simulations start. Indeed, a larger vari¬ 
ance in the reionization histories on these scales is seen by 
comparing c orrespondingly sized sub-volumes of larger-scal e 
simulations (llliev et al.ll2006l : llliev et al.l2014| : lGnediij2014l l. 
However, this variance is caused primarily by the large range 
of environments sampled by the sub volumes, each of which 
may have a different mean density lllliev et al.ll2014l ~). The 
top right panel of Figure [2] shows further that there is a sig¬ 
nificant dependence of the reionization history on resolution. 
At higher resolution, reionization ends at higher redshifts. 
This mostly reflects the increased SFRs at increased resolu¬ 
tion that our simulations show during the early and middle 
phases of reionization. 

In our simulations, reionization completes at higher red- 
shifts than suggested by constraints on the neutral frac¬ 
tion from observations of Lyman-a emitt ing galaxies and 
high- redshift quasars (e.g., Figure 5 in [Robertson et al.l 
l2013ll . However, interpretations of these observations are 
dominated by l arge systematic uncer t ainties (for a dis¬ 
cussion see, e. g_^ iRobertson et al.l l2013l: IJe nsen et ^ l2013l : 
iDiikstral l2014l : iMesinger et al.l 20 14 Tavlor fc Lid d l2014l l . 
At 2 = 6, the simulated mean neutral fractions are in 
reasonable agreement with observatio nal constraints from 
the Lyman-a forest ([Fan et a l.| |2006l), which are thought 
to be relatively robust (e.g., F^net^ 200^ but see, e.g., 


to be relatively robust I e.g.. Itan et al.M2UUa: but see, e.g., 

McGreer. Mesinger. fc Fanll201ll:lGnedinll2014:[Becker et alJ 


2014h . Nevertheless, it is important to keep in mind that 
our simulations may underestimate the duration and over¬ 
estimate the redshift of reionization because the limited size 
of th e simulated volume may b i as th e reionization histories 
(e.g., [iliev et ^ [ 2 OI 4 I : [GnedirJ [2014[ l. and because of our 
choice of a large sub-resolution escape fraction = 1-0. 

A useful check for consistency between the total SFR 
densities and the reionization histories is provided by com¬ 
parison with the critical SFR density Psfr nee ded to keep 
th e IGM ionized (lMadnUjJaardtj_^£^,eea|l99^ ; updated as 
in [Pawlik. Schave. fc van Scherpenzeell [ 200 ^ 1 . 


PSFR « 0.013 M 0 yr ^ Mpc ® 

/ GiGM ) f fesc 


1 + 2 


( 6 ) 
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Figure 1. Images of the gas overdensity, the ionized fraction, the temperature, and the neutral fraction (left to right, top to bottom) 
in the reference simulation L25N512 at z ^ 9. Each image is an average in a central slice of linear extent 25 h~^ cMpc and thickness 
2.5 h~^ cMpc cut perpendicularly to one of the axes of the simulation box. Reionization is already well underway (the volume-averaged 
ionized fraction is ~ 0.4) and has heated the gas in ionized regions to temperatures ~ 2 X 10^ K. The gas near galaxies may reach 
temperatures up to ~ 10^ K as it is heated by SNe and structure formation shocks. Our spatially adaptive RT simulations reveal a 
network of self-shielded neutral clumps and filaments inside the ionized regions. 


where /esc is the fraction of ionizing photons 
that, on averag e, escape galaxies to io nize the 


IGM(e.g., |Rjazpumoy J^SpnjrrieM^ 


Gnedin^ Kravtsov, fc Chen 


2008 


Paardekooper. Khochfar, fc Dalla VecchiJ 


Wise fc CenI 


200e; 


200C: 


2013; 


Yaiima, Choi, fc Nagamind l201lh . and ~we have assumed 


a Chabrier IMF, which implies an ionizing emissivity 
per unit stellar mass larger, and therefore a critical SFR 
density smaller, by a factor of about 1.7 compared to that 
implied by a Salpeter IMF (see Section 12.31) . An IGM 


clumping factor Cigm = 5 is a conservative estimate from 
our sim ulations (see Figure El), and consistent wi th earlier 
works llPawlik, Schave. fc van Scherpenzee I I 2 OO 9 I : see the 
summary of results in Finlator et al. I 2 OI 2 I ). 


In this case, the simulated star-forming galaxies can 
keep the gas at z < 7 ionized provided /esc ^0.3. Whether 
such a relatively large escape fraction is physically plau¬ 
sible is a matter of active research (see, e.g ., the recent 
overview in lBenson, Venkatesan. fc Shul]||2013l) . Computing 
the escape fraction from our simulations would require ad- 
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Figure 2. Effect of box size and resolution on the star formation and reionization histories. Curves of different colors correspond to 
different simulations as indicated in the legend. Top left: comoving SFR density. Dashed and solid curves show the total SFR density 
and the SFR density including only the contribution from galaxies with SFRs greater than 0.3 Mq yr“^, the current detection limit 
(assuming a Chabrier IMF, corresponding to a limiting UV magnitude Mab see Figure ! ^. The dotted curve shows the critical SFR 

Eg. [^ , assuming CjciM/fesc = 25. Observation al data is taken fr om|Bouwens et al.| |2pl4a|), iFinkelstein et al.l ll2 014[l. [Bouwens et al.l 
201^) . lDuncan et alj||2014l . their SEP models). | Oesch et al] ll2013l) . rOesch et^dTI^^l^ . lCoe et al.lll2013^ . and lEllis et al.l ll2013l . revised 


down by a factor of two following lOesch et al.|[2013h . as indicated in the legend, and divided by 1.7 if necessary to convert from the 
Salpeter IMF to the Chabrier IMF used here. When necessary, we have corr ected the published SF R densities upwards to account for 
contribution from galaxies down to the current detection threshold, following iBouwens et all l|2014al'l . Top right: volume-weighted mean 
neutral hydrogen fraction. The horizontal dotte d line marks a neutral fraction = 0.5 to guide the eye. The (model - depen dent) 
constraints on the neu t ral fra ction are taken f romlSch roeder^ Mesinger^ Haiman| l|2013l . QSO damping wingsh fOuch^_et__aI.I |2^^, Ly-a 
clustering), lOta et ^ 1 1 I 2008 I. Ly-a emitters), iDiikstra, Mesinger, &: Wvithj ll201ll . Ly -a emission fraction) and iFar^^alTll20^ . Ly-a 
forest), following the discussion of Figure 5 in fRobertsoi^t al.l ll2013lv Bottom left: mean (solid) and median (dashed) volume -weighted 
hydrogen photoionization rate. The triangle with error bars marks the observational constraint from [Wvithe &: BoltonI ll201ll') . Bottom 
right: electron scattering optical depth towards reionization. The horizontal dotted line shows the latest estimate from observations by 
the Planck satellite, and the grey band indicates the associated 1-sigma error interval. 
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ditional RT computations beyond the current work. Obser¬ 
vations currently p robe SFRs > 0.3 Mp yr~^ (assumiuR a 
Chab rier IMF; e.g.. 


I 2 OI 4 I : iDuncan et al 


Bcn^ens et al.ll^l4al . iFinkelstein et all 


2014 Figure El). Under the above as¬ 


sumptions, at z = 6, the observed population of galaxies 
alone is capable of keeping the Universe ionized, lending 
support to scenarios in which reionization is driven by star¬ 
form i ng galaxies (e.g., Pawlik. Schave ^ fc van Scherpenzeell 


I 2 OO 9 I : lArihelstein et ahl 20121 : ICai et al.il2014h . 

The simulations predict a strong increase in the volume- 
averaged mean and median hydrogen photoionization rates 
as the simulation volumes are reionized (bottom left panel 
of Figure O. This signals a rapid build-up of an ioniz¬ 
ing background, which sets the equilibrium neutral hydro¬ 
gen abundance in the r eionized IGM le.g.. I^nedinl |200^ : 
lAubert fc TevssiCTl2010f) . The median photoionization rate 
is in good agreement with constraints at z = 6 from ob¬ 
serva tions of, e.g., the Lyman-a f orest llB olton fc Haehneltl 
|2007^ and quasar near-zones le.g.. IWvithe fc BoltorJ 201lll . 
However, the mean photoionization rate exceeds the me¬ 
dian by more than an order of magnitude. The excess 
in the mean photoionization rate with respect to the ob¬ 
servational constraints is similar to that seen in previous 
simulations of reionizatio n (e.g., lAubert fc TevssiCT I 2 OIOI : 
iFinlator. Dave. Ozelllioilh . 

The optical depth towards reionization, Treion = 
riecdt, where rie is the number density of free elec¬ 
trons and Zstart ^ 30 is the redshift at which reionization 
begins, is an impo rtant integral constraint on the reioniza- 
Alvarez et al. I I 2 OO 6 I : IShull fc VenkatesarJ 


I 2 OO 8 I : lAhn et ^|2012| ~). For the computation of the optical 
depth, we assume a hydrogen mass fraction X = 0.75, and 
that the fraction of ionized hydrogen is equal to that in our 
simulations, which assume X = 1. Furthermore, we make 
the standard assumptions that hydrogen is fully ionized at 
z < 6, that helium is neutral (singly ionized) whenever hy¬ 
drogen is neutral (ionized) at z > 3, and that helium is 
doub ly-ionized at z < 3 le.g.. Illiev. Scannapieco. fc Shapii^ 

[200i). 

The bottom right panel of Figure 0 shows that 
our simulations yield optical depths consistent with the 
most recent constraints derived from observations of the 
cosmic microwave background, Treion = 0.066 ± 0.016 
llPlanck Collaboration et al.ll201^ . Our simulations lack the 
resolution and the physics to follow the formation of stars 
in minihaloes, i.e., haloes with virial temperatures below 
10"* K. Ionization by these star s is expected to increase 
the optical depth by < 0.03 (e.g.. IShull fc Venkatesanll200§ : 
lAhn et ^l2012l : IWise et al.ll2014h . If these stars evolve into 
accreting black holes, the optical depth may be further in- 
creased due to the pre-io nization of the IGM by X-rays (e.g., 
iRicotti fc OstrikCTl20o4l . If the additional contribution to 
ionization from minihaloes and X-ray sources is indeed that 
high, and if we were to add this contribution to our estimates 
of the optical depth, then this would yield values in slight 
excess of the Planck observational estimates. However, con¬ 
sistency with observations may still be achieved by lowering 
the ISM escape fraction of ionizing photons, which is a free 
parameter in our simulations. In this case, reionization in 
our simulations would occur later and yield a lower optical 
depth (see Figure [H. 

Current cosmological simulations of reionization often 



Figure 3. Effect of stellar feedback on the SFR history in the 
reference simulation L25N512. Shown are SFR densities in the 
simulations without feedback (noFB; red dashed), with only SN 
feedback (noPH; red solid), with only photoheating (noSN; blue 
dashed), and with both SN feedback and photoheating (L25N512; 
blue solid). SNe strongly reduce the cosmic SFR. On the other 
hand, photoheating has a comparatively small impact on the cos¬ 
mic SFR, and this impact is nearly negligible in the absence of 
SNe. 


have difficulties in simultaneously matching observational 
constraints on the neutral fraction, the photoioniza¬ 
tion r ate and the optical depth (e.g., lAubert fc Teysshm 
l201fll : IFinlator. Davi. Ozell l2nilh . Figure [2] shows that 
our simulations are no exception. Many works have 
argued that these difficulties may signal a significant 
role by physical processes that are currently not under¬ 
stood or unknown, and that are currently only poorly 
modelled or lacking altogether in cosmological simula¬ 
tions. Such processes include, among others, an evolving 


IGM 

2 OI 2 I: 


escape fraction (e.g., Kuhlen fc Faucher-Gigufera 


AlrargZjFinlgigrj^.^gMj 


2 OI 2 I : 


Mitra. Ferrara, fc GhoudhurvI l2013l : Ferrara fc LoebI l2013h 


and the absorption of photons by Lyman limit s ystems 
fe.g.. [^lton fc Haehneltll2ni3l : iMeknger et ahlli^nidl : which 
our reference simulation begins to resolve). Some works 
have successfully matched the observati onal constraints b y 
a careful calibration of parameters fe.g.. ICiardi et aDl2012l ). 
In Section ED we confirm that calibrating the ISM escape 
fraction of ionizing photons helps to improve the match with 
observational constraints. However, it is also important to 
acknowledge that current observational constraints on the 
neutral fraction at z > 6.0 are highly model-depen dent and 
too weak to single o ut reionization scenarios fe.g.. [piikstral 
|2014 iGnedinI I20l4l . Investigating these issues is beyond 
the aims of the current work. Here we will focus on how 
feedback from photoheating and SNe impacts the formation 
of galaxies and the reionization of the IGM in physically 
motivated scenarios of reionization. 
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Figure 4. Evolution of the UV luminosity function in the ref¬ 
erence simulation L25N512 between 2 = 6 — 9 (solid curves) . 
Symbols show observational constraints from [Smit^t_ah|_|201^ 
at 2 = 6 (filled squares, corrected for dust) and Duncanetaf 


a open squares, correct ed for dust), fron i iBo uwens et~.. 
, open diamonds) andjFinlglstein etjrlJ ||2014| . filled dia¬ 


mo nds) at 2 = 8 and f rom lOescli^r'ahn^Ol.'j . open triangles) 
and lMcLeod et al] ||2014| . filled triangles, their estimate assuming 
density evolution) at 2 ^ 9. The comparison at 2 fa 7 is included 
in Figure[5]for clarity. We have omitted any upper limits reported 
in these works. Dashed and dotted curves are the corresponding 
Schechter fits, as indicated in the legend. The simulated UV lu¬ 
minosity function is in good agreement with the observational 
constraints across the entire range of redshifts. 


3.2 Stellar feedback 

In this section we investigate how stellar ionizing radiation 
and SNe impact the formation and evolution of galaxies and 
the reionization of the IGM. We will compare the simula¬ 
tions listed in Table [T] which include both photoheating and 
SNe, with simulations in which photoheating and SNe are 
disabled but which are otherwise identical. We use the suf¬ 
fix noFB to distinguish the latter simulations from the for¬ 
mer. We further compare with simulations that include SNe 
but not photoheating, distinguished by the suffix noPH, and 
with simulations which include photoheating but not SNe, 
distinguished by the suffix noSN. We focus on our reference 
simulation L25N512 to illustrate the physics at play, and we 
assess the impact of box size and resolution on our conclu¬ 
sions. 


3.2.1 Impact on the cosmic SFR history 

Figure[3]quantifies the impact of stellar feedback on the evo¬ 
lution of the total cosmic SFR density in comparisons of our 
reference simulation L25N512 with simulations that differ 
only in the inclusion of the stellar feedback processes. The 
comparison with simulation L25N512-noFB, in which both 


photoionization heating and SNe were turned off, reveals the 
strong suppression of star formation by stellar feedback. The 
comparison with simulations L25N512-noPH, in which only 
photoheating is turned off, and L25N512-noSN, in which 
only SNe are turned off, shows that the suppression of the 
SFR due to SNe is much stronger than the suppression of 
the SFR due to photoheating. 

SNe have a significant impact already at early times, 
2 < 10, long before the IGM in the simulation is substan¬ 
tially ionized. On the other hand, the impact of photoheat¬ 
ing becomes noticeable only at 2 < 8, when the average ion¬ 
ized fraction approaches unity. The belated impact of pho¬ 
toheating on the SFR suggests a stronger role of suppression 
of star formation due to the illumination of galaxies by the 
ionizing background than due to the exposure to local ion¬ 
izing radiation from stars inside the galaxies. However, the 
impact of photoheating on the cosmic SFR remains compar¬ 
atively small also after reionization. 

Figure [3] also shows that photoheating suppresses star 
formation more strongly in the simulation with SNe than in 
the one without SNe (compare the difference between the 
two solid curves with the difference between the two dashed 
curves). On the other hand, SNe suppress star formation 
more strongly in the simulation with photoheating than in 
that without (compare the difference between the two blue 
curves with the difference between the two red curves). Pho¬ 
toheating and SNe thus amplify each other in suppressing 
the formation of stars. These results are in goo d agreement 
with our earlier work in IPawlik fc Schavd (l2009l ). in which we 
first reported this effect in simulations in which the gas was 
in photoionization equilibrium with a uniform ionizing back¬ 
ground instantaneously turned on at 2 = 9, and in which SN 
feedback was implemented by kicking gas particles in winds 
(see also Finlator. Dave. Ozell I2OIII : iHambrick et al.l 1201 ll : 
iHopkins et al.l 20131 ). 


3.2.2 Impact on the UV luminosity function 


Another important observable of galaxy formation is the 
UV luminosity function, i.e., the volume density of galaxies 
per unit logarithmic UV luminosity. It quantifies the con¬ 
tributions to the cosmic SFR density from galaxies with 
different SFRs, and therefore helps us to understand the 
processes that shape the cosmic SFR density discussed 
above. Figure |4] shows the UV luminosity function at 
2 « 6 — 9 in our reference simulation L25N512. We con¬ 
verted SFRs and UV luminosities using SF R/( Mf;i yr~^) = 
0.6 X 1 . 25 X 10~^^Luv/( erg s“^ Hz“^) (e.g.. lBouwens et aJl 
l2014a: iKennicu'^ Il99a ). which yields the SFR function 
(e.g., Smit et al. 20121 ). The factor 0 .6 converts the origi¬ 


nal relation, which is based on a Salpeter IMF, to a re¬ 
lation assuming the Ghabrier IMF used here (see Sec¬ 
tion [2]3]). Luminosities are related to AB magnitudes using 
—2.5logjQ[Luv/(47rlO^ pc^)] — 48.6. Our reference simula¬ 
tion L25N512 yields a SFR function in good agreement with 
the observations above the limiting SFR of > 0.3 M© yr“^ 
(assuming a Ghabrier IMF) currently probed by observa¬ 
tions. Small difference exist at 2 > 9, where observational 
estimates are characterized by large uncertainties, and at 
2 < 7, especially at SFRs > 2 x 10 M© yr“^, where dust 
corrections may be significant. 

In the left-hand panel of Figure [5] we investigate how 
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Figure 5. Comp arison of simulated and observed SFR functions a,t z ^ 7. Squar es show the observed d ust-corrected SFR function from 
ISmit et alj ||2012|) . which is based on the dust-uncorrected SFR functions from [Bouwens et alj ll201lh . We have divided the observed 
data, which assumes a Salpeter IMF, by 1.7 to enable comparisons with our simulations, which assume a Cha brier IMF. The grey d otted 
curve shows the corresponding Schechter fit. Diamonds show the observed dust-corrected SFR function from lPuncan et al.l J 2 OI 4 I . their 
SED models). Left: effect of stellar feedback in the reference simulation L25N512. The top panel shows the ratio of the SFR functions 
in the bottom panel and the SFR function in the absence of stellar feedback extracted from simulation L25N512-noFB. Feedback from 
SNe is critical to matching the observed SFR function. The inefficiency of SNe in suppressing star formation at low SFRs is primarily a 
numerical artefact of the limited resolution. Photoheating affects the SFR function at SFRs ^0.2 Mq yr~^, but this is largely masked 
by the impact of SNe, which is strong at SFRs > 0.05 Mq yr“^. Photoheating thus does not imprint a noticeable characteristic scale in 
the UV luminosity function. Right: dependence on box size and resolution. At higher resolution, the explosion of massive stars as SNe 
suppresses star formation more strongly, although the effect is relatively small and therefore all curves are close to each other and to 
observational estimates. The turnover in the SFR functions at low SFRs is primarily a result of the limited resolution and the lack of 
low-temperature gas physics. 


feedback impacts the UV luminosity function, where we fo¬ 
cus for simplicity on a single characteristic redshift z ^ 7. 
The comparison of the reference simulation with simulation 
L25N512-noFB, in which both photoheating and SNe are 
turned off, shows that suppression of star formation by stel¬ 
lar feedback is a critical element of simulations that attempt 
to match the observed SFR function. In our reference simu¬ 
lation, reionization heating and SNe suppress star formation 
at nearly complementary scales. The SFR function in simu¬ 
lation L25N512-noSN, in which only SNe are turned off and 
which serves to demonstrate the isolated impact of photo¬ 
heating, approaches that in simulation L25N512-noFB at 
SFRs > 0.2 Mq yr“^. On the other hand, the SFR function 
in simulation L25N512-noPH, in which only photoheating is 
turned off and which serves to demonstrate the isolated im¬ 
pact of SNe, approaches that in simulation L25N512-noFB 
at SFRs < 0.05 Mq 

The inefficiency of photoheating at high SFRs is ex¬ 
pected if SFRs correlate positively with halo mass, as is the 
case in our simnlations (Figure |6}, since re ionization affects 
most l y the gas in low-mass haloes (e.g., Barkajm^fcj^bl 
I 1999 I : IPilkstra et al.|[2004l : lOkamoto. Gao, fc Theun^ 2008ll . 


However, the decreased efficiency of SN feedback at fow 
SFRs is likefy a numericaf artefact of the fimited reso¬ 
lution. Indeed, we show in Appendix m that in galaxies 
with low SFRs, star formation is more strongly suppressed 


by SNe at higher resolution. Note that our simulations 
also likely underestimate the local impact of photoheating 
by stars, especially in the lowest-mass star-forming galax¬ 
ies (iRahmati et al.l l2013lJl , thus likely raising the relative 
importance of photoionization by the cosmological back¬ 
ground. 


Dependence on resolution and box size 

The right-hand panel of Figure [5] shows how our results de¬ 
pend on resolution and box size. While our reference sim¬ 
ulation yields SFR functions in good agreement with the 
observations, simulations at lower and higher resolution sig¬ 
nificantly overpredict and underpredict the observed SFR 
functions in the range of SFRs in which SN feedback is ef¬ 
ficient. This can be understood by noting that the rate at 
which the SN heated g as cools and loses energy dep ends 
on resolution (Eq. 17 in iDalla Vecchia fc Scha^l2012lf . At 
higher resolution, the SN heated gas can maintain a high ra¬ 
tio of the local radiative cooling time and the sound-crossing 
time across a resolution element, a requirement for efficient 
feedback, out to higher densities, increasing the abifity of 
the SN heated gas to prevent the formation of new stars. 
Because at higher resolution, the gas fraction in low-mass 
haloes is more strongly reduced by stellar feedback, also the 
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hierarchically assembling larger-mass haloes grow systemat¬ 
ically mo re baryon-deficient as reso lution is increased (see 
also, e.g.. iFinlator. Dave. dzelll201lh . 

The differences between the simulated and observed 
SFR functions in the high-res and low-res simulations may 
be reduced by increasing or lowering the fraction /sn of 
the energy each SN injects. This would be reasonable since 
our simulations lack the resolution and the physics that is 
needed to accurately model the radiative losses in the ISM. 
This strategy would enable one to investigate numerical con¬ 
vergence in simulations that reproduce observational con- 
straints on star formation indep endent of resolution (e.g., 
lGnedinll2014l : ISchave et al.l[201^ . The current set of sim¬ 
ulations, on the other hand, highlights the dependence on 
resolution of the SN energy fraction required to match ob¬ 
servational constraints on star formation. This dependence 
demonstrates clearly that current cosmological simulations 
have limited power in predicting the properties of star¬ 
forming galaxies from first principles. Note however, that 
the difference in the simulated and observed SFR func¬ 
tions is relatively small, of the order of the difference in¬ 
troduced by uncertainties in the conversion between UV lu¬ 
mino si ties and SFRs (see discussions in, e.g., Munoz&Loed 


20jT];_jBo YjmnKojnl^ Bullock, fc Garrison-Kimmell 12014 


Bouwens et al.l l2014al f . Therefore, the agreement between 
simulated and observed SFR functions is still very good. 

Finally, we investigate the dependence on box size. 
The reference simulation overpredicts the abundance of the 
most intense star-forming galaxies in comparison with ob¬ 
servations. This is likely due to the finite size of the box, 
which biases estimates of the abundance of the most mas¬ 
sive and therefore rarest galaxies. The comparison with the 
small and large box simulations L12N256 and L50N512 
supports box size as a limiting factor in accurate deter¬ 
minations of the SFR function at the highest SFRs. Ad¬ 
ditional processes that may impact models of the SFR 
function at high SFRs include feedback from massive 
black holes (e.g., Di Matteo. Springel. fc HernauistI l2005l : 


iBooth fc Schavgl2009ll , which we have ignored here, and ob - 
scuration by dust (e.g., ICen fc Kimm|[2014l : ICai et aLll2014h . 
While we compare to observations of dust-corrected SFRs, 
there are large systematic uncertainties in the amount of 
dust as well as its pro perties, especially at the high redshift s 
of interest here fe.g.. fsmit et al.ll2012l : IPuncan et al.ll2014h . 
It is also important to keep in mind that observational es- 
timat es of the SFR function are subject to cosmic variance 
(e.g., iFinkelstein et al.|l2012f) . 


show different SFRs (e.g.. [^rkana fc LoebI[ 2 OO 0 II . because 
of hierarchical assembly that leaves also relatively massive 
haloes deficient in baryons as they assemble from low-mass 
photoevaporated progenit ors (e.g., iBarkana fc LoebI I 2 OOC 1 I : 
IFinlator. Dave. C)zelll201ll l. and because of the averaging in 
volumes i n which reionization m ay proceed in a patchy man¬ 
ner (e.g., IBarkana fc Loebll2006f) . 

In our simulation L25N512-noSN, which does not in¬ 
clude SN feedback and therefore lets us isolate the im¬ 
pact of radiative heating, photoheating indeed causes a pro¬ 
nounced change in the slope of the SFR function at SFRs 
< 0.2 Mq yr“^. However, in simulation L25N512 that in¬ 
cludes both photoheating and SN feedback, the change in 
the slope of the UV luminosity function due to photoheat¬ 
ing is strongly masked by the suppression of star forma¬ 
tion by SNe at SFRs > 0.05 Mq yr“^. In the presence of 
SNe, the change in the slope of the SFR function occurs 
at lower SFRs, near the scale at which the lack of low- 
temperature physics prevents gas cooling and star forma¬ 
tion in our simulations. This result is robust to changes in 
resolution explored here. We thus conclude that if SN feed¬ 
back is as efficient as suggested by our simulations, detect¬ 
ing the signature of reionization heating in the SFR func¬ 
tion will be challenging. However, it is important to keep in 
mind that our simulations ignore that stars may also form in 
minihaloes cooling through molecular hydrogen transitions, 
which would be more susceptible to feedback from photo¬ 
heating. 

The steep rise of the SFR function in the presence of 
radiative heating and SNe at low SFRs down to the SFRs 
that our simulations do not accuratel y resolve is consistent 
with that seen in previous works (e.g., Finlator. Dave. Ozell 


S revious works (e.g., Fi 
2 OI 2 I : iGai et al.ll2014ll . and with ex trapo¬ 
lations of the observed SFR function (e.g., Smit_et_^ 20li 


Bmyens et al.|[2014al : IFinkelstein et al.ll2014l : Duncan et alJ 


2014) . Because the mass scale below which photoheating af¬ 
fects the SFRs of galaxies is close to the mass scale below 
which efficient gas cooling requires the presence of molec¬ 
ular hydrogen or metals, quantifying the precise impact of 
reionization on the UV luminosity function will require more 
detailed, higher-resolution simulations of galaxy formation 
that include the relevant gas physics. While this is chal¬ 
lenging because it extends the relevant dynamic range to 
still smaller scales, impressive first ste ps towards such simu¬ 
lations have already been ma de (e.g.. iHaseeawa fc Semehiil 
I 2 OI 3 I : iGnedin fc Kauro^l20l4 . 


The UV luminosity function as a probe of reionization 

The reduction of the SFRs of low-mass galaxies due to 
photoheating may affect the slope of the SFR function at 
low SFRs. Several (semi-) analytical works have proposed 
to search for this type of signature of reionization in the 
SFR function (e.g.. lBarkana fc Loebll^OOl : iMashian fc LoebI 
|2014 . Because the mass scale below which reionization sup¬ 
presses star formation depends on the timing of reionization, 
locating t his signature may constrain the reionizatio n his¬ 
tory (e.g. . IBarkana fc Loebliooil : iMufioz fc Loebll201ll l. Un¬ 
fortunately, the imprint of photoheating in the SFR function 
may be quite weak and therefore difficult to detect, because 
of intrinsic scatter that causes galaxies of the same mass to 


3.2.3 Impact on galaxy properties 

Inspecting the properties of individual galaxies helps gain 
insight into the physical origin of the shape of the SFR 
function discussed above. Figure |6] shows the median SFRs 
(top row) and the median baryonic mass fractions /bar = 
(M* -I- Mgas)/Mvir (bottom row) of galaxies in our simu¬ 
lations. The left-hand panel in each row shows the impact 
of feedback in the reference simulation at z « 7, and the 
right-hand panel shows the evolution with redshift. In the 
top right panel, we have divided the SFRs by the virial mass 
of the haloes hosting the galaxies to improve the clarity of 
the presentation. Our reference simulation predicts that the 
faintest galaxies accessible by current observations reside in 
haloes with masses Mvir ~ 10^° Mq (intersection of solid 
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Figure 6. Median SFRs (top) and median baryonic mass fractions /bar (bottom) of galaxies of v irial masses Mvir and associated circular 
velocities = (GMvirAvir)^/^ = 17(M/10® M0)l/3[(1 + z)/10\^/'^ km (e.g., Eq. 3.11 in iLoebllioToll . In each row, the left-hand 
panel shows the impact of stellar feedback in the reference simulation L25N512 at a typical redshift 2 7 and the right-hand panel shows 

the dependence on redshift in that simulation. In the right-hand panels, to isolate the impact of photoheating, we include results from 
simulation L25N512-noSN, in which SN feedback was turned off (dashed curves). In the top right panel, we have divided the SFRs by 
the virial masses of the haloes hosting the galaxies to improve the clarity of the presentation. The dashed vertical lines mark the masses 
of haloes with virial temperatures T^ii. = 10'* K. The dotted vertical lines mark the mass of 100 DM particles. The dotted diagonal line 
in the top right panel marks a SFR of 0.3 Mq yr“*, which is the limiting SFR accessible in current observations (assuming a Chabrier 
IMF; see Figure[4]|. The dashed horizontal lines in the bottom row panels mark the cosmic baryon fraction f!i,/(f^b + f^DM)' The strong 
suppression of star formation by photoheating and by the explosion of stars as SNe is accompanied by a strong reduction in the baryon 
fraction. The dependence on resolution is discussed in Figure EH 


and diagonal dotted curves in the top middle panel), con¬ 
sistent with observational estimates based on matching the 
shapes of the DM halo m ass function and t he observed UV 
luminosity function (e.g.. iTrenti et aLll2010l ). 

The top left panel of Figure [6] shows that in our ref¬ 
erence simulation, the median SFR at 2 « 7 drops sharply 
below Mvir ~ 2 X 10® M©. The comparison with simulation 


L25N512-noSN in which SNe were turned off, shows that this 
scale is set primarily by feedback from reionization. In the 
absence of photoheating, star formation continues in haloes 
with lower masses. However, even without feedback, galax¬ 
ies in haloes with masses < 5 x 10® M© do not form stars 
since gas cooling and condensation is prevented by the finite 
resolution and the lack of low-temperature coolants, such as 
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Figure 7. Effect of SNe on the evolution of the volume-weighted 
mean ionized fraction and its dependence on resolution. The 
solid curves show the reionization histories in simulations includ¬ 
ing feedback from ionizing radiation and SNe. The dashed curves 
show the reionization histories in the corresponding simulations 
in which SN feedback was turned off. The inclusion of SNe leads 
to at most a small delay in reionization, despite the strong sup¬ 
pression by SNe of the cosmic SFR (FigurelSj. 


molecular hydrogen and metals. SN feedback suppresses star 
formation strongly across nearly the entire range of galaxy 
masses. It only becomes inefficient at the lowest masses at 
which our finite resolution impacts the conversion of gas into 
stars, and at the highest masses, at which our simulations 
are also subject to finite box size effects. 

Stellar feedback reduces star formation primarily be¬ 
cause it reduces the amount of dense gas, both by expelling 
gas and by limiting the rates at which gas is accreted. In 
our reference simulation, reionization strongly reduces the 
gas fractions at halo masses < 2 x 10® M©. Reionization 
also reduces, though less strongly, the gas fractions in up 
to ~ 10 times more massive galaxies, which assemble in 
mergers of lower-mass baryon-defici ent galaxies and accrete 
gas from the reionized IGM (e.g., Barkana fc Loelj l2000l : 
iFinlator. Dave. Ozell l201ll : iMunoz fc LoebI feed¬ 

back most strongly impacts the gas fractions at masses 

> 10 ® Mq in our reference simulation, again demonstrating 
that reionization and SN feedback act mostly complemen¬ 
tary at the reference resolution. We show in Appendix 
that in our high-res simulation, both reionization and SN 
feedback reduce the baryon fractions more strongly, and SN 
feedback extends its impact to haloes with masses as low as 

> 2 X 10® M0. 

The scale at which photoheating becomes effective in 
suppressing star formation and reducing the baryon frac¬ 
tions in our reference simulation evolves significantly only 
after 2 « 8, which coincides with the completion of reioniza¬ 
tion. A similar effect is seen in the evolution of the baryon 
fractions. This late impact of photoheating suggests a larger 
role of illumination by the external ionizing background 


than by localized ionizing sources internal to the galaxies, in 
agreement with the discussion of the belated impact of ra¬ 
diative heating on the cosmic SFR density above. However, 
we caution that in our simulations the impact of local ion¬ 
izing sources may be underestimated due to the finite reso¬ 
lution. At least before reionization, radiative feedback from 
local sources is expected to dominate over feedback from 
external illumination. This is indeed seen in simulations of 
the first stars and the start of reionization, which focus on 
smaller volumes and therefore can afford a higher resolution 
and also accurately capture the relevant low-temperature 
physics, primar ily the cooling and chemistry of molecular 
hydrogen (e.g., Rjcotjd_j^^jriec^^ 20051: Wis£_^^Te| 20081 : 

iM 


IPawlik. Milosavlievic. fc Brornni 20131: 1 Jeon et al.l 2015ll. 


3 . 2.4 Impact on reionization 


Both SNe and radiative heating may have a strong im¬ 
pact on reionization. We have seen above that SN feed¬ 
back reduces the SFRs and therefore the ionizing luminosi¬ 
ties of galaxies, which makes it more difficult for the galax¬ 
ies to reionize the Universe. However, SNe may also open 
low-density channels in the ISM, e.g. by expelling gas in 
winds, thro ugh which ionizing photons may escape more 
easily (e.g.-lDove . Shull, fc Ferraril2000l: i Fuiita et S]l2003l : 

enl 


IWise fc CerJ 20091 : Paardekooper ^ al.l 20111 '). This makes it 

easier for galaxies to reionize the Universe. The net impact 
of SNe on reionization is the result of the interplay of these 
two processes, and is therefore difficult to predict. 

We have also seen above that radiative heating may help 
SNe to suppress star formation, which impedes reionization. 
On the other hand, radiative heating raises the Jeans mass 
in the IGM, and this reduces the IGM clumping factor and 
therefore the rate at which the IGM recombines. If recom¬ 
binations consume a signihcant number of the ionizing pho¬ 
tons that escape into the IGM, then radiative heating will 
help keeping the ionized gas ionized and thereby facilitate 
reionization (e .g., Pjwylik^_Schaye^_jzvaj_ScheQ3emee]||2^ 
Finlator et al.l I 2 OI 2 I : ISo et al.l l2014l : ISobacchi fc Mesingerl 


2014| . 


Figure [3 shows that SN feedback has only a small im¬ 
pact on the timing of reionization. In the reference and low- 
res simulations, the reionization histories with and without 
SN feedback are nearly identical, and in the high-res simula¬ 
tion, SN feedback delays reionization by A 2 < 0.5. That the 
net impact of SNe on reionization is small suggests that the 
strong reduction in the SFRs due to SNe - in the high-res 
simulation by a factor ~ 10 by the end of reionization (Fig¬ 
ure 0 - is partially compensated by an increase, also due to 
SNe, in the escape fraction of ionizing radiation. A plausible 
physical mechanism by which this might be achieved, is the 
strong reduction of the baryon fractions by SNe seen in our 
simulations (Figure 0. However, the lack of a substantial 
impact on the reionization history could also be explained if 
reionization is driven mostly by the lowest-mass galaxies in 
which SN feedback is less effective in suppressing star for¬ 
mation due to the limited resolution. RT computations of 
the escape of ionizing photons into the IGM would help to 
identify which of these mechanisms is dominant, but this is 
beyond the scope of the current work. 

The left-hand panel of Figure [8] shows that 
the reduction in the IGM recombination rate due 
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Figure 8. Evolution of the IGM clumping factor CioOi which parametrizes the average recombination rate in gas with overdensity ^ 100. 
Left: effect of stellar feedback in the reference simulation L25N512. Photoheating strongly reduces the clumping factor because it increases 
the Jeans mass in the reionized gas, providing a positive feedback on reionization. SN feedback increases the clumping factor by ~ 15% 
as it moves gas from galaxie s to the IGM. Right: dependence on box si ze and resolution. The dashed curve shows the clumping factor 
from simulation r9L6N256 in iPawlik. Schave. fc van Scherpenzeell ll2009h . in which the gas was heated by a uniform ionizing background 
turned on at 2 ^ 9. The earlier result is in good agreement with the clumping factor in simulation L12N512, which employs a similar 
resolution and in which the IGM, on average, is reionized at a similar redshift. The small increase with respect to the earlier work is 
mostly due to the inclusion of SNe. 


to radiative heating is strong. We have followed 
IPawlik. Schave. fc van Scherpenzeell (l2009l 'l and have 
equated the IGM clumping factor, Cigm, to Cioo, which 
parametrizes the average recombination rate of gas with 
overdensities ^ 100. This enables us to separate re¬ 

combinations in the IGM from recombinations insid e 
galaxies (see also iMiralda-Escude. Haehnelt. fc Reesll2000ll . 
The latter are typically parametrized by the ionizing 
escape fraction. Alternative definitions of the clumping 
factor are sometimes employed, e.g., by applying ad- 
ditio nal selection criteria to identify the ionized IGM 
fe.g.. [K ohler. Gnedin, fc Hamiltonl 1200 ^: IShull et al.l l2012l : 


iFinlator et al. 20121 : Kaurov fc GnedirJ 2014 ). However, for 
typical reionization scenarios, these definitions, which are 
all designed to achieve the same objective, i.e., locating the 
ionized gas in the IGM, ge nerally yield very si milar results 
(see, e.g., the discussion in IFinlator et al.]l2012l l. 

Our simulations likely underestimate the clumping fac¬ 
tor before reionization, since the Jeans mass in the unheated 
IGM is unresolved (e.g., lEmberson. Thomas, fc Alvared 
l2013h . Thus, the positive radiative feedback from photoheat¬ 
ing on rei onization is also underestimated. Not e howe ver, as 
shown in IPawlik. Schave. fc van Scherpenzeell (l2009l ). that 
the clumping factor at 2 « 6 is insensitive to the redshift at 
which reionizat ion occurs for reionization at 2 > 8 . As a l- 
ready found in IPawlik. Schave. fc van Scherpenzeell ll2009ll . 
SNe move gas from galaxies to the IGM, and this leads to a 
slight increase in the clumping factor (see also lFinlator et ^ 

l2012lb 

The right-hand panel of Figure [8] shows that the 
clumping factor derived from the high-res simulation 


L12N512 is consistent with that in simula tion r9L6N256 of 
IPawlik. Schave. fc van Scherpenzeell (120091 '). which had sim¬ 
ilar resolution but assumed that the IGM is heated by 
a uniform UV background turned on instantaneously at 
2 = 9. This agreement likely results because in our RT sim¬ 
ulations, reionization occurs sufficiently rapidly such that 
differences in the time at which individual regions inside 
the simulated volume are reionized are small and simi¬ 
lar to the time it takes for the IGM to respond dynam¬ 
ically to the increase in the Jeans mass in the simula¬ 
tion in which the gas is instantaneously exposed to a uni¬ 
form ionizing background. The slight increase in the clump¬ 
ing factor by ~ 15% with respect to the earlier work is 
mostly du e to the inclusion of SN feedback, a lthoug h, as 
shown in IPawlik, Schave. fc van Scherpenzeell (l2009ll . the 
differences in cosmological parameters also contribute. For 
reionization occurring at redshifts 2 > 8, the IGM clump¬ 
ing factor Cioo at 2 « 6 is numerically converged and in- 
sensitive to a further increase in resolut ion and box size 
dPawlik. Schave. fc van Scherpenzeell [200^ ~) . 


4 DISCUSSION 

In this section we briefly discuss how changes in physical 
parameters affect the outcome of our simulations and also 
mention some of the main physical processes our simulations 
have ignored. 
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Figure 9. Impact of variations in the physical parameters on the volume-weighted mean neutral fraction (left), the hydrogen photoion¬ 
ization rate (middle), and the median SFRs of galaxies at z ~ 7 (right, normalized by the virial masses) in simulation L12N256. The 
vertical dashed line marks the mass corresponding to a virial temperature 10^ K, and the vertical dotted line marks the mass of 100 DM 
particles. Data points are as in Figure [2] Reducing the sub-resolution escape fraction to 0.6, 0.5, and 0.4 (blue, green, and light green 
curves) implies that reionization occurs later, which reduces the impact of photoheating on the SFRs. Increasing the energy injected in 
each photoionization by a factor 4 (red curves) leads to larger gas temperatures, which reduces recombination rates and increases the 
scale below which photoheating suppresses SF. 


4.1 Variations of physical parameters 


Figure [9] shows how variations in the physical parameters 
impact the reionization history and the ability of galaxies to 
form stars. Towards this end we have repeated the small-box 
simulation L12N256, which has the same resolution as the 
reference simulation L25N512 but allows for a computation¬ 
ally more efficient exploration of the parameter space, with 
different values for some of the parameters. 

One of the main parameters of reionization simulations 
is the escape fraction of the unresolved ISM. Decreas¬ 

ing this fraction from 1.0 to 0.6, 0.5 and 0.4 causes reion¬ 
ization to occur at lower redshifts, as expected, at which 
it extends over longer tim es, in good agreement with sim- 
ilar parameter studi e s by IPetkova fc Springell (1201 ll) and 
iHasegawa &: SemelinI (l2013f ). Since photoheating is delayed, 
low-mass haloes can continue to form stars more efficiently 
down to lower redshifts. However, the associated increase 
in the ionizing emissivity is smaller than the reduction due 
to the decrease in the sub-resolution escape fraction, and 
so the impact on reionization due to the change in the 
minim um mass of haloes is small (e.g., IPetkova fc Springell 
l201lh . In the current simulations, adopting an escape frac¬ 
tion ~ 0.4 yields excellent agreement with observa¬ 

tions of the evolution of the ionized fraction and the hydro¬ 
gen photoionization rate at 2 = 6 (but observational con¬ 
straints are weak; see the discussion in Section [3.111 . 


The amount of energy injected in each photoioniza¬ 
tion depends on the spectrum of the radiation sources 
and requires multi-frequency RT simulations for an ac- 


Abel & Haehneld 

1999bl: iMaselli. Ciardi. & Kanekail l2009l: 

Pawlik & Schavd 

2011h. Here we treat this energv as 

a oarameter (e.e.. Petkova & Sorineell |201l|). Increas- 


ing it implies a slight acceleration in reionization, and 
a slight decrease in the mean neutral fraction after 
reionization. This is caused by the increase in the 
gas temperatures implied by the higher photoheating 


rates, which in turn dec reases the rate at which hy¬ 
drogen recombines (e.g., |SHuyellLJ)aJl,_J^Jkmagia| |2()04|; 


PayHik, Schave. fc van Scherpenzeell l2009l : iFinlator et al.l 


2 OI 2 I) . Finally, the increased gas temperature increases the 


negative impact of photoheating on the efficiency of low- 
mass galaxies to form stars, raising the mass scale below 
which star formation is s trongly suppressed (see also, e.g., 
IPetkova fc Springell201ll ). 

Finally, we have carried out a preliminary compari¬ 
son with results from a new set of simulations of reioniza¬ 
tion similar to those presented here, which feature increased 
physical realism and span a wider range of box sizes and 
resolutions. These simulations will be discussed elsewhere 
in more detail and include, among others, helium chem¬ 
istry and cooling/heating, and feedback and metal enrich¬ 
ment from AGB stars, core-collapse and Type la SNe, and 
metallicity-dependent population synthesis. The treatment 
of helium was made feasible by replacing the explicit chem¬ 
istry solver described in Sec. m. with the implici t solver 
described in lPawlik, Milosavlievic. fc BrommI (l2013l ). which 
is faster. Moreover, the simulations are designed, by cali¬ 
brating the sub-resolution escape fraction and the SN energy 
fraction, to match the observed UV luminosity functions and 
to exhibit similar reionization histories independent of res¬ 
olution. A simulation employing the same resolution as our 
reference simulation here and adopting a sub-resolution es¬ 
cape fraction of ~ O-Sj yields reionization and SFR 

histories similar to those in the current simulation L12N256 
adopting an escape fraction ~ 0-6- 


4.2 Limitations 

Our simulations ignored a number of potentially relevant 
physical processes. Most importantly, perhaps, our simu¬ 
lations ignored the chemistry of and radiative cooling by 
molecular hydrogen, effectively assuming a soft UV back¬ 
ground that prevents the build-up of hydrogen molecules. 
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This approximation fails at the earliest stages of reion¬ 
ization, where it artificially prevent s the formation of 
stars i nside low-m a ss mi ni haloes (e.g., JWisej_rurkj_&_^be| 
2008; Greif et ^ 20081: JPawlik, Milosavlievic. fc BrommI 


2013 : [Muratov et al. 2013li . This early population of stars 

may preionize the IGM and provide a significant feed- 
back on subsequent star formation and reionization (e.g., 
iRicotti fc Qstril^l2004^ : lAhn et ^l2012l i . Our simulations 
have also ignored the enrichment with metals that accom¬ 
panies the explosion of stars as SNe. Metal-enrichment 
affec ts the rates at which gas cools and forms star s 
(e.g.. [jappsen et '^120091 : 1 Wiersma. Schave. fc Smit'hll2009lf . 
which may be especially important in low-mass minihaloes 
in which radiative cooling by atomic hydrogen is sup¬ 
pressed. Finally, we have only followed the radiative feed¬ 
back from the relatively soft ionizing radiation emitted 
by metal-enriched stars. Other sources of radiation, such 
as zero-metallicity stars (Pop III stars) or X-ray emit¬ 
ting black holes, may be an important source of ionization 
and feedback during reionization (e .g., Madau_gt-aT 2004 
Ricotti. Ostriker. fc GnedinI l2005l: lAlvarez. Wise, fc Abell 


20091 : HaimanI 201 ll : Jeon et al. l2014ab 


5 SUMMARY 

We have carried out a suite of cosmological radiation- 
hydrodynamical simulations of galaxy formation during 
reionization. The reference simulation was run in a box of 
size 25 h~^ cMpc and contained 512® dark matter and 512® 
baryonic particles, thus resolving atomically cooling haloes 
with at least > 10 dark matter particles. Simulations using 
both larger and smaller boxes and higher and lower resolu¬ 
tion allowed us to investigate the numerical convergence of 
our results. Simulations in which either SNe or photoioniza¬ 
tion heating or both are turned off, enabled us to isolate and 
investigate the impact of feedback from star formation. Ion¬ 
izing photons were transported using accurate and spatially 
adaptive RT, tracking the growth of ionized regions and the 
build-up of an ionizing background at the native high reso¬ 
lution at which the hydrodynamics was carried out. 

Current cosmological simulations lack both the resolu¬ 
tion and the physics to provide an ab initio description of 
the structure of the interstellar gas and the rate at which the 
gas cools. This necessitates the use of physically motivated 
but resolution-dependent parameters to control the energy 
that each SN injects and the fraction of ionizing photons 
that escape into the IGM. SFRs and reionization histories 
are sensitive to these parameters and this impedes the use of 
cosmological simulations in predicting these quantities from 
first principles. On the other hand, one may exploit this 
sensitivity and choose parameters for which simulated SFRs 
and reionization histories are consistent with current obser¬ 
vational constraints, and investigate the implications of such 
observationally supported models of galaxy formation. Here 
we have focused on the role of feedback from SNe and photo¬ 
heating, two processes that critically shape galaxy formation 
and reionization. 

Our reference simulation yields SFR densities and a 
UV luminosity function in excellent agreement with obser¬ 
vational constraints and completes reionization by 2 « 8. 
Increasing the resolution leads to a strong increase in the 


cosmic SFR at high redshifts as it facilitates the conden¬ 
sation of gas into low-mass galaxies. It leads to a mild de¬ 
crease in the cosmic SFR at late times, when star formation 
is strongly regulated by stellar feedback. As a consequence, 
near the end of reionization, our high-res simulation yields a 
slightly smaller normalization of the UV luminosity function 
and our low-res simulation yields a slightly larger normal¬ 
ization of the UV luminosity function than our reference 
simulation and observations. Because of the initial increase 
in the cosmic SFR, increasing the resolution also increases 
the redshift at which the IGM is reionized. Increasing the 
size of the simulation box above 12.5 cMpc has only a 
minor impact on the SFR and reionization histories. 

Photoheating reduces the baryon fractions and sup¬ 
presses star formation primarily in haloes with masses below 
< 2 X 10® Mq. SNe, on the other hand, reduce the baryon 
fractions and suppress star formation primarily in haloes 
more massive than > 10® M©. Therefore, the currently ob¬ 
servable cosmic SFR is more strongly suppressed by SNe 
than by photoheating, and SN feedback alone is sufficient to 
match observational constraints on the UV luminosity func¬ 
tion. The inefficiency of SNe in the lowest mass galaxies is 
primarily a consequence of the lack of low-temperature gas 
physics and the limited resolution. Nevertheless, the feed¬ 
back from SNe is sufficiently strong to mask the impact of 
photoheating on the abundance of low-mass atomically cool¬ 
ing star-forming galaxies. We thus do not find a noticeable 
signature imprinted by reionization heating on the UV lumi¬ 
nosity function, although we note that the resolution of our 
simulations is insufficient to model star-forming minihaloes 
that would be more strongly affected by photoheating. 

Despite the relatively small impact on the cosmic SFR, 
photoheating is an important feedback process. First, pho¬ 
toheating amplifies the ability of SNe to suppress star for¬ 
mation. This amplification is nonlinear and mutual, demon¬ 
strating the need to simultaneously account for both feed¬ 
back processes. Second, photoheating smooths out gas den¬ 
sity fluctuations in the IGM and thereby strongly reduces 
the IGM recombination rate. This makes it easier to keep 
the gas ionized, which facilitates reionization. In contrast, 
the net impact of SNe on reionization is small. SNe strongly 
suppress SFRs and slightly increase the IGM recombination 
rate as gas is moved from the galaxies to the IGM. How¬ 
ever, this leads only to a small delay in the timing of reion¬ 
ization, possibly because SNe create additional low-density 
channels in the ISM through which ionizing photons can 
escape, which increases the escape fraction of ionizing radi¬ 
ation, or because reionization is driven by the lowest mass 
galaxies in which SN feedback is inefficient in our simula¬ 
tions. 

In summary, our work demonstrates that both photo¬ 
heating by the stellar radiation that reionizes the Universe 
and the explosion of massive stars as SNe may have had a 
strong impact on structure formation and reionization in the 
first billion years. 
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APPENDIX A: DEPENDENCE ON 
RESOLUTION 

Here we extend the discussion of the galaxy properties in 
our reference simulation in Section 13.2.31 (Figure [6l) with 
a brief investigation of the dependence on resolution. Fig- 
ure l All shows that in our high-res simulation, in the absence 
of SN feedback, reionization suppresses star formation only 
at masses below ~ 5 x 10® Mq. This is significantly less 
than at the reference resolution, even though the IGM in 
the two simulations is reionized at similar times. The de¬ 
crease in the suppression scale is caused by the increase in 
the SFRs that, in the absence of feedback, accompanies an 
increase in resolution, and for which photoheating does not 
entirely compensate in our simulations. However, the high- 
res and reference simulations agree closely on the character¬ 
istic scale ~ 10® Mq below which the median SFR is strongly 
suppressed by the combined action of radiative and SN feed¬ 
back. In the low-res simulation, star formation is suppressed 
in haloes as massive as 5 x 10® Mq, primarily by the limited 
resolution. 
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Figure Al. Impact of resolution on the median SFRs (left) and the median baryonic mass fractions /bar (right) of galaxies of virial 
masses Mvir at 2 : 5^ 7. The meaning of the dashed and dotted lines is as in Figure [6] 
















